Correlations in a BEC Collision: First-Principles Quantum Dynamics with 150 000 Atoms 
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The quantum dynamics of colliding Bose-Einstein condensates with 150 000 atoms are simulated directly 
from the Hamiltonian using the stochastic positive-P method. Two-body correlations between the scattered 
atoms and their velocity distribution are found for experimentally accessible parameters. Hanbury Brown-Twiss 
or thermal-like correlations are seen for copropagating atoms, while number correlations for counterpropagating 
atoms are even stronger than thermal correlations at short times. The coherent phase grains grow in size as the 
collision progresses with the onset of growth coinciding with the beginning of stimulated scattering. The method 
is versatile and usable for a range of cold atom systems. 
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The prediction of many-body quantum dynamics is a 
long term goal of investigation in a variety of scientific 
fields ranging from physics to chemistry, biology and com- 
putation theory. It is a pivotal problem for interacting 
systems, but challenging because of the complexity of a 
full description of a quantum system, in which the num- 
ber of basis states grows exponentially with the number 
of particles. Experiments with Bose-Einstein condensates 
of ultra-cold atoms give excellent examples of phenom- 
ena that are not well described by standard approximations 
such as the Gross-Pitaevskii (GP) equation. This equation 
treats the macroscopically occupied wavefunction, but ne- 
glects atomic correlations and fluctuations!]]] which are es- 
pecially prominent in strongly interacting or dimensionally 
reduced gases, and in condensate collisions. In the latter 
case, the GP equation fails because the scattering initially 
occurs spontaneously into unoccupied modes, which are ig- 
nored by a macroscopic wavefunction approach. Later, the 
scattering becomes Bose-enhanced, and a coherent, non- 
perturbative treatment of the scattered modes is essential. 

Treatments of BEC collisions have included a slowly- 
varying envelope approximation (SVEA) which esti- 
mates the scattering cross-section^, perturbation theoryyl 
01, and the semi-classical truncated Wigner methodj5|]. 
The last method is nonperturbative, works well in one 
dimension[6], and appears to treat both the initial sponta- 
neous scattering and the later Bose enhancement. However, 
we will show that it gives strongly incorrect results in 3D 
at large momentum cutoff, because the equations of motion 
are truncated. Hence, there is a strong incentive to develop 
a quantitative, first-principles method for these cases. 

This Letter also has a broader focus than just BEC. While 
path-integral Monte Carlo methods are now very successful 
for calculating equilibrium properties, quantum dynamics 
is not amenable to these techniques because of the very 
rapid dephasing between different paths[7]. Phase-space 
distribution methods (such as the Glauber-Sudarshan|8], 
positive-P|0], stochastic wavefunction[10], gauge-P ll lln 
do not suffer from this problem, and yet the scaling is still is 
only linear in the system size. They have been applied suc- 



cessfully to cold atom quantum dynamics in increasingly 
large systems, including simulation of evaporative cool- 
ing to form a BEC|12], spin squeezing and formation of 
two-component BECs lU3ll . correlation dynamics in a uni- 
form gas[ 14], the quantum evolution of Avogadro's number 
of interacting atoms|15], the dynamics of atoms in a ID 
trap lfioll . and molecular down-conversion lfl^l . 

Here we demonstrate the maturity and ready-to-use na- 
ture of the original positive P method for truly macroscopic 
systems. We simulate an average of 150 000 atoms, requir- 
ing M = 1.08 x 10 6 momentum modes. Since each of 
M modes can have up to about N atoms , the full Hilbert 
space contains at least D w M N w 10 1 000 000 orthogonal 
quantum states (or D w io 200000 if fixed total atom num- 
ber is assumed). This is one of the largest Hilbert spaces 
ever treated in a first principles quantum dynamical simula- 
tion — made possible by probabilistic sampling rather than 
brute-force diagonalization. 

The use of such a first principles, yet stochastic, simu- 
lation confers several advantages in comparison with ap- 
proximate methods. Firstly, all uncertainty in the results 
is confined to random statistical fluctuations, with no sys- 
tematic bias. This uncertainty can be reduced by averag- 
ing over more stochastic realizations, and even more im- 
portantly, can be reliably estimated from their spread. Sec- 
ondly, these methods lead to relatively simple equations of 
motion, which are easily adapted to realistic modeling of 
trap potentials and local losses. 

We consider the collision of two pure 23 Na BECs, with 
a similar design to a recent experiment at MIT lfl7ll . A 
1.5 x 10 5 atom condensate is prepared in a cigar-shaped 
magnetic trap with frequencies 20 Hz axially and 80 Hz ra- 
dially. A brief Bragg laser pulse coherently imparts a ve- 
locity of 2vq = 19.64 mm/s to half of the atoms, much 
greater than the sound velocity of 3.1 mm/s. At this point 
the trap is turned off so that the wavepackets collide freely. 
In a center-of-mass frame, atoms are scattered preferentially 
into a spherical shell in momentum space with mean veloci- 
ties v s as vq. As the density of atoms in this shell builds up, 
Bose-enhancement of scattering into it is expected to begin. 
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Of particular interest are the distribution of scattered atom 
velocities, and correlations between those atoms, which 
were recently shown to be experimentally measurable lfl8ll . 

In present BEC experiments, the system can be described 
to a high accuracy by the local interaction Hamiltonian lfl9ll : 
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The operator W(x) creates a bosonic atom at position x = 
(x, y,z) and obeys commutation relations ^'(y)] = 

S ^ (x — y), with 5 a delta-function tempered by a momen- 
tum cutoff fc max . The coupling constant g depends on the 
s-wave scattering length a (2.75nm in the case of 23 Na), and 
for fc max Cl/d one finds that g = 4irH 2 a/m. 

To calculate time-evolution, we employ the positive P 
representational [Till because it preserves the full quan- 
tum dynamics. This approach utilizes the completeness 
of the coherent-state basis[8]. The density matrix p is 
expanded as a positive distribution P over off-diagonal 
coherent-state projectors [9], thus preserving quantum cor- 



relations: p = J P{aJ)d 2M ad 2M p \a)((3*\/[(/3*\a}}. 
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operators in a volume V. The coherent state \a) = 
®£ \ a k) * s a J omt eigenstate of each aj, with complex 
eigenvalue a gig]. When used to expand the master equa- 
tion ihdp/dt = [H,p] (3, this leads to a Fokker-Planck 
or diffusion equation in the probability P, which is equiva- 
lent to solving an ensemble of stochastic equations for the 
sampled variables a and (3. The equations are simplified on 
discrete Fourier transforming to a conjugate spatial lattice 



as (2) 
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Here V 2 as is the discretized analogue of V 2 a(x) for a 
field, and and £g are real Gaussian noises, independent 
at each time step (of length At) and lattice point, with stan- 
dard deviations 1/VAVAt, where AV = V/M . 

There is an equivalence between statistical averages 
of moments of aj: and and corresponding normally- 
ordered expectation values of operators a^, at. As the 
number of trajectories, S, grows towards oo, the correspon- 
dence becomes exact. These stochastic equations are just 
the mean-field GP equations in a doubled phase-space, plus 
noise terms. Remarkably, these modifications incorporate 
all effects beyond the GP equation, provided certain phase- 
space boundary conditions are met ! 1 ll 12011 . 

Uncertainty in the observables is estimated by binning 
trajectories, then calculating the observable predictions 
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Figure 1: Evolution of velocity distribution. The distribution 
p(v x ,v y ) has been integrated over one transverse dimension z. 
Color from blue to red indicates increasing density (its range varies 
between panels). S — 2048 trajectories. 



from each bin, and using the central limit theorem to esti- 
mate the standard deviation in the final mean of bin means. 
Lattice spacings At and Ax are chosen by reducing them 
until no further change is seen. In the figures, results 
are presented in terms of velocity space, v = hk/m, the 
Fourier transformed field ^f(v), and the velocity space den- 
sity p(v) = (&(v)^(v)). 

Following earlier procedures !^ , we discretize onto a 
M = 432 x 50 x 50 lattice with fc x , max = 1.4 x 10 7 /m 
and fcy^maic = 6.2 x 10 6 /m. We begin the simulation in 
the center-of-mass frame at the moment the lasers and trap 
are turned off (t = 0). The fc max and lattice size are cho- 
sen large enough to encompass all relevant phenomena but 
small enough that the spacing (7r/fc max ) is much larger than 
a. The initial wavefunction is modelled as the GP solu- 
tion of the trapped t < condensate, but modulated with 
a factor [e jfcQX + e ~ ifc Q x ] j which imparts initial veloc- 
ities v x = ±vq = ±hkQ/m in the x direction. For com- 
putational reasons, the mean number of atoms in the sys- 
tem is 1.5 x 10 5 here, compared to « 3 x 10 7 in the MIT 
experiment! 17]. As in other recent treatments(5|], we ig- 
nore thermal atoms and initial quantum depletion, for sim- 
plicity. For our parameters, 10% thermal component will 
occur at ss 0.38T c lll], giving a ss 1% quantum depletion of 
the ground state in the center of the cloudily. These small 



corrections can be included in the initial state[23]. 

Figure Q] shows the formation of the scattered atom shell. 
Careful inspection shows that the mean scattered atom 
speed v s is less than the wavepacket speed |vq|, as noted in 
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Figure 2: Correlations between scattered atoms. All are be- 
tween scattered atoms at a maximum density point in the shell 
(with velocity v = (0, -9.37, 0) mm/s relative to the COM) and 
those with a shifted velocity v, where: a: v — (0,v y ,0), b: 
v = (v x , 0, 0), and c: v = (0, 0, v z ). To reduce statistical noise, 
averages g± — J gr'™' (yo + 8v, v ± 5v) d 2 8v, over a volume 
Vq in velocity space are plotted l28ll . Triple lines are la errors. 
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Figure 3: Evolution of correlations between scattered atoms. 

Parameters and ?7o as in Fig. [2] b: The Full-width half-maximum 
(FWHM) width of \g± | in velocity space. 

J3|. We also see weak scattering between two atoms from a 
wavepacket at ±vq to one atom at ±3vq and one at Tvq. 

Ranged two-body correlations give insight into typi- 
cal small-scale behaviour during a single experimental 
run. The first-order correlation function g^(vi,V2) = 
(vi)^ (v^)) I \J p{v\) p(v2), describes coherence be- 
tween particles with velocity vi and v%. The second- 
order (number) correlation function g 1 - 2 ' (vi, V2) = 
(^^Vl) 1 i ji (v2)^(vl) 1 i | (v2)}/p(vl) p{v2) gives the average 
shape and size of "lumps" in the velocity distribution. 

The dynamics of the correlations among scattered atoms 
are shown in Figs.[2]and[3] Locally the atoms are thermally 
bunched with g^ 2 '(v, v) w 2 in a "Hanbury Brown-Twiss" 
manner (Fig. |2}. This behaviour has been confirmed qual- 
itatively in a similar recent He* experiment [24]. The lo- 
cal region over which coherence is strong, dubbed a "phase 
grain" by Norrie et al\5], is described by \g^\. It closely 



matches the condensate wavepackets' p(v) in size, and is 
wider than g^ by w V2 (Fig. |2j». We find that the orienta- 
tion of these phase grains is constant throughout the whole 
spherical shell. Interestingly, after t as 200 /is, the phase 
grains expand significantly in the radial direction (relative 
to COM) (Figs. [2k and[3})). This onset of growth coincides 
with the beginning of Bose stimulated scattering (see below 
and Fig.@h, circles). 

Atoms with velocities v and — v on opposite sides of the 
spherical shell are not coherent (|g < - 1 -'| ~ 0), but are corre- 
lated in number (Fig. [30. Initially, correlations are extreme: 
g( 2 '(k, —k) ^> 2. This is analogous to a two-mode mixed 
state p 1 1, 1) (1, 1 1 + (1 -p) |0, 0) (0, 0| with a small probabil- 
ity p of single atoms in both modes and otherwise vacuum. 
There g^> = 1/p, At longer times, g^ is seen to decay 
in Fig. [3^, although it is still much greater than the thermal 
value of two for t > 200/is when the phase grain contains 
several atoms. To measure short time velocity correlations, 
one might try to preserve them by suddenly switching off 
the atomic interactions using a Feshbach resonance during 
the collision. After expansion, they would develop into po- 
sition correlations fisll . 

Some previous correlation estimates are in qualitative 
agreement: For longer times, g^ (v, v) — 2, as well as 
g^(v,—v) f» 2 and g^(v,—v) ~ were predictedyD. 
Truncated Wigner calculations^ saw the presence of phase 
grains, but their orientation or dynamics were not studied. 
High initial correlations may have not been seen due to the 
known poor signal-to-noise ratio in that method. 

The scattering rate (Fig.Hk, circles) goes through two dis- 
tinct phases: The spontaneous regime of constant scattering 
into almost empty modes is seen for 30 ps < t < 200 ps, 
followed by the stimulated (Bose-enhanced) regime for 
times t > 200 ps, where there is a decided increase in scat- 
tering rate despite a lessening overlap between the collid- 
ing wavepackets. We interpret this transition as the onset 
of Bose enhancement of scattering into the spherical shell 
around |u| w v s = 9.37 mm/s. As a rough check, it should 
begin when the number of particles in a locally coherent re- 
gion ("phase grain") approaches one. Using the widths of 
from Fig.[5J>, and the calculated density at |uj = v s one 
finds ks 0.9 atoms per phase grain at 200 ps. 

A comparison to approximate methods used previously 
is instructive. Fig. H shows our total predicted scatter- 
ing rate and the distribution of axial (x) velocities, com- 
pared to the approximate truncated Wigner method. The 
accuracy of that method is very poor with our parameters, 
even surprisingly so. It adds a halo of false particles (de- 
tail in Fig H)>) out to about ±2vq, while at higher veloc- 
ities unphysical negative densities are obtained. Since it 
is a hidden-variable theory it must introduce half a virtual 
particle per mode in the initial conditions to model vac- 
uum fluctuations, but it does not distinguish them from 
the "real particles". Then, virtual particles at velocity v 
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accessible with this approach. 
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Figure 4: Comparison of exact and approximate methods. 

Panel a: Rate of scattering out of coherent wavepackets. Obtained 
by counting atoms outside of spheroidal regions covering the co- 
herent wavepackets centered at ±vq and ±3vq. Panel b: Distri- 
bution of axial (x) velocity of scattered atoms at t = 657/is; In 
both panels: S = 2048 (Pos. P), S = 672 (Wigner). 

are scattered by the condensates at rs ±vq in the process 
v & ± vq — > v' & (v± vq — v'). As a result, modes at 
high velocities become depleted compared to the physical 
vacuum, while the extracted virtual particles accumulate at 
lower velocities and take on the appearance of a real density, 
as was also discussed previously J25|]. 

For any single momentum mode this effect is small, but it 
becomes very significant when a large number of modes are 
calculated. A relatively higher momentum cutoff will in- 
crease the error, as the fraction of virtual particles increases 
(or vice-versaf!]). This indicates a generic ultra-violet di- 
vergence of the error with the truncated Wigner method. 

The main limitation of the positive-P method is the 
growth of sampling uncertainty with time. It eventually 
reaches a size where it is no longer practical to produce 
enough trajectories for useful precision. In our case useful 
results are obtained for t < 660 /is. This useful time range 
depends on several factors, with coarser lattices, weaker in- 
teractions, or smaller density all extending it lll4il . Signif- 
icant extensions appear achievable by tailoring appropriate 



stochastic gauges lll 11 12CX1 or basis sets to particular systems. 

In conclusion, we have simulated the quantum dynamics 
of macroscopic interacting Bose gases from first-principles, 
obtaining momentum space densities and ranged correla- 
tion functions for atoms scattered during the collision of two 
BECs. Previous approximate calculations were partly veri- 
fied, while a variety of new phenomena are also predicted, 
including the growth of phase grains in the radial momen- 
tum direction, and strong correlations at short times be- 
tween scattered atom pairs. The truncated Wigner method 
was confirmed strongly incorrect in regimes where the num- 
ber of condensed atoms per lattice site is less than one. 

This demonstrates that phase-space methods are a tool 
that is ready-to-use for first principles calculations for ex- 
perimentally realizable systems. Similar calculations ap- 
pear feasible for a broad range of cold atom systems (in- 
cluding fermions lferjIO . A range of other phenomena that are 
difficult to describe quantitatively with approximate meth- 
ods (e.g. macroscopic EPR and entanglement ll27l0 may be 
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